Evidence for out-of-equilibrium crystal nucleation in suspensions of oppositely 

charged colloids 
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We report a numerical study of the rate of crystal nucleation in a binary suspension of oppositely 
charged colloids. Two different crystal structures compete in the thermodynamic conditions under 
study. We find that the crystal phase that nucleates is metastable and, more surprisingly, its 
nucleation free energy barrier is not the lowest one. This implies that, during nucleation, there is 
insufficient time for sub-critical nuclei to relax to their lowest free-energy structure. Such behavior is 
in direct contradiction with the common assumption that the phase that crystallizes most readily is 
the one with the lowest free-energy barrier for nucleation. The phenomenon that we describe should 
be relevant for crystallization experiments where competing solid structures are not connected by 
an easy transformation. 
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Liquids often must be cooled substantially below the 
freezing temperature before spontaneous crystallization 
occurs in the bulk. The reason is that the system has 
to overcome a free-energy barrier when moving from 
the metastable liquid to the stable solid phase. When 
the two phases are separated by a high free-energy bar- 
rier, spontaneous fluctuations that would result in the 
formation of the stable phase are unlikely and there- 
fore rare. Most fluctuations will result in the formation 
of ephemeral "sub-critical" crystal nuclei that redissolve 
spontaneously. Only occasionally a crystal nucleus will 
form that exceeds the critical size needed for spontaneous 
subsequent growth [l[ . The crystal nucleation rate is de- 
fined as the number of post-critical clusters that form 
per unit time in a unit volume. In classical nucleation 
theory (CNT) , it is assumed that sub-critical clusters are 
in quasi-equilibrium with the parent phase This as- 
sumption is reasonable if the time it takes to establish an 
equilibrium distribution of sub-critical clusters is short 
compared to the time needed to nucleate a crystal. If 
the nucleation rate is low, the steady-state distribution 
of sub-critical clusters of size n is (nearly) proportional 
to exp(— j3AG(n)), where AG(n) is the free energy asso- 
ciated with the formation of a crystalline cluster of size 
n in the metastable liquid. The CNT expression for the 
nucleation rate per unit volume is 



(1) 



where k is a kinetic prefactor and AG cr n is the height 
of the nucleation barrier. It was already pointed out 
by Ostwald 0] that often, during crystal nucleation, a 
solid phase forms that is not the thermodynamically most 



stable one. Stranski and Totomanow [|[ have rational- 
ized this observation in the language of CNT by sug- 
gesting that the phase that nucleates is the one sepa- 
rated from the parent phase by the lowest free-energy 
barrier - and this need not be the most stable solid 
phase. Implicit in this explanation is the assumption that 
the kinetic prefactor k is similar for different nucleation 
routes, and that hence the relative nucleation rates are 
exclusively determined by the heights of the nucleation 
barriers. Unfortunately, the assumptions underlying the 
Stranski- Totomanow (ST) rule cannot easily be tested in 
experiments. Here we present simulations where we com- 
pute independently the rate of crystal nucleation and the 
height of the free energy barriers separating a metastable 
liquid from two more stable solid phases. 

In order to study rare events such as liquid-solid nu- 
cleation by simulation, one has to resort to special simu- 
lation techniques, precisely because a typical nucleation 
event does not occur within the time scale of a conven- 
tional simulation. The only alternative is to use very 
large system sizes Q and long simulation times [Bj]. But 
even then the metastable system has to be prepared in 
a state deeply supersaturated before spontaneous nucle- 
ation can be observed . 

Here, we use the the Forward-Flux-Sampling (FFS) 
method of Allen et al. 0, @] to compute the rate of crys- 
tal nucleation. This method was designed to study rare 
events both in and out of equilibrium. It can be used 
under conditions where brute-force simulations become 
impractical. FFS has been used to calculate the rate of 
crystal nucleation in molten sal ts M and the nucleation 
rate of an Ising model in pores [1 01 ] . 

To compute the free-energy barriers for crystal nucle- 
ation, we use umbrella sampling 11]. This method has 
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been used before to compute the free-energy barriers for 
the nucleation of crystals 12, IH and liquids [3]. The 
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umbrella-sampling approach determines the variation of 
the free energy of the system with a reaction coordinate 
that measures the progress of the transformation from 
the liquid to the crystalline phase. It should be noted 
that, whereas the nucleation rate is an observable quan- 
tity, the height of the free-energy barrier for crystal nucle- 
ation may depend somewhat on the choice of the reaction 
coordinate. 

Making use of the information obtained using both 
methods we will show below that crystal nucleation in 
a mixture of oppositely charged colloids is incompatible 
with the ST conjecture. 

In our simulations, we studied a 1:1 binary mixture of 
monodisperse, oppositely charged colloids. The screened 
Coulomb interaction between two colloids of diameter a 
and charge Ze is approximated by a Yukawa potential: 
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where the sign is positive for equally charged and neg- 
ative for oppositely charged colloids, Xb — e 2 /e s kBT 
is the Bjerrum length (e s is the dielectric constant of 
the solvent) and k = \f&nXBPsait is the inverse Debye 
screening length (p sa it is the number density of added 
salt). A hard core prevents colloids from overlapping. 
The total energy of the system is the sum of the pair 
interactions. The cut-off radius, r c , is 3.5cr. We define 
the reduced temperature T* = (1 + na /2) 2 a / Z 2 Xb as the 
inverse of the contact energy, and the reduced pressure 
as p* = pT*a 3 /ksT. The phase diagram of this poten- 
tial for ka = 6 reproduces the solid structures that are 
found experimentally in mixtures of oppositely charged 
colloids [15[ . In this system, two solid phases can coexist 
with the fluid. At high temperatures, the liquid phase 
coexists with a substitutionally disordered face-centered 
cubic colloidal crystal (disordered- fee). At low tempera- 
tures the stable solid at coexistence has CsCl structure, 
where the charges are ordered on a bec lattice (Fig. [T|) . 

This system is a suitable candidate to test the Stranski- 
Totomanow conjecture as two distinct solid phases may 
form during crystal nucleation. In contrast to systems 
that have been studied earlier [3, [13] these two solids 
are not connected by an "easy" (e.g. martensitic) trans- 
formation. Besides, we study nucleation close to the co- 
existence temperature between both solids. 

Both for the FFS calculations and for the calculation 
of the free-energy barrier separating liquid and solid, we 
need a reaction coordinate that measures the progress 
of the nucleation process. In the present study, we use 
n, the number of particles in the largest solid cluster, 
as a reaction coordinate [l8j |. Our reaction coordinate 
distinguishes liquid from solid but is not sensitive to the 
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FIG. 1: Phase diagram of the system under study in the T*, 
packing fraction (r) = -k<t 3 N/6V) plane [l5j]. L+S stands 
for liquid-solid coexistence. The circle indicates the state 
point of the metastable liquid (T*=l, p*=15, ?7=0.526).Tfie 
disordered-fcc-fluid-CsCl triple point temperature is T* = 
1.07. 



structure of the crystal lattice. In fact, under the thermo- 
dynamic conditions used in our study (T* = 1, p* = 15), 
every particle in either the disordered fee phase or the 
CsCl-like solid is identified as crystalline. On the con- 
trary, in the metastable liquid phase, less than five out 
of 1000 particles were identified as crystalline. 

The FFS technique expresses the nucleation rate, 4>i- s , 
as the product of two factors [19( : 
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where <f>i—j is the rate at which spontaneous fluctuations 
lead to the formation of a small crystallite consisting of J 
particles, whilst Pj- S denotes the probability that such 
a cluster will grow to form a bulk solid, rather than re- 
dissolve. In what follows, we will ignore the effect of 
hydrodynamic interactions and estimate <pi—j using a ki- 
netic Monte Carlo (MC) algorithm [20j with a maximum 
displacement of 0.01 a [Is |21|. In the limit of small trial 
displacements, the MC algorithm approaches Brownian 
Dynamics, but MC has the added advantage that we can 
easily work in the NpT ensemble. The probability Pj- S 
is computed as a product of probabilities: 
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where Pj-k is the probability that a trajectory that 
starts with a cluster of size J, will grow to size K rather 
than redissolve. This probability can be estimated by 
starting a number of trajectories from a cluster of size 
J and dividing the number of those that arrive at K by 
the total number of trials. The successful trajectories 
provide starting configurations for the next step, namely, 
the calculation of the probability that cluster K will grow 
to size L, rather than redissolve. The FFS method only 
works if the dynamics of the system is not fully determin- 
istic. In the present case, different kinetic MC trajecto- 
ries (with maximum displacement 0.04 cr) were generated 
from the same configuration by changing the seed of the 
random number generator. We stress that the "reaction 
coordinate" in the FFS scheme is only used to measure 
the progress of the crystal growth -it does not favor one 
crystal structure over another 0, Q . 
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In this work, we have studied the crystallization of the 
metastable liquid phase in a system of 1000 particles at 
T* = 1 andp*=15 (p*/p* oex ~ 1.7). The packing fraction 
of the liquid at the coexistence pressure, p* ocx = 8.8±0.1, 
is rj = na 3 N/6V = (0.471 ± 0.005). At T* = 1 and 
p* = 15, the packing fraction of the metastable liquid is 
(0.526 ± 0.005). 

At p* = 15, no spontaneous nucleation is observed even 
after 3- 10 6 MC cycles (a cycle consists of a trial move per 
particle and a volume move). However, the nucleation 
rate can be computed at p* — 15 using the FFS method 
([3]). We find that the probability that a crystalline clus- 
ter of 5 particles will continue to solidify is 10 _28±2 . The 
rate at which spontaneous fluctuations in the metastable 
liquid result in the formation of crystalline clusters of 5 
particles is equal to 10 _4±1 Docr -5 (where Do denotes the 
diffusion coefficient at infinite dilution). Therefore, the 
nucleation rate is estimated to be 10 _32±3 Z?oc r 5 - We 
observe the growth of a rather compact solid cluster in 
the metastable liquid. When the crystalline cluster has 
reached a size of n = 120± 15, it has a 50% probability of 
redissolving: this is our operational definition of the crit- 
ical nucleus size. To identify the crystal structure of the 
solid nucleus, we analyse the radial distribution function 
of the particles that belong to the cluster. This provides 
a convenient way to distinguish disordered fee and CsCl 
structures. Figure [2] a shows the comparison of g(r) for a 
cluster of 80 particles with that of the bulk solid phases. 
From the figure, it is evident that the arrangement of 
the particles in the growing solid cluster is fcc-like even 
though the stable solid phase at T* = 1 and p* = 15 is 
CsCl. 

The formation of crystal nuclei of a metastable solid 
phase can be interpreted as a manifestation of the Ost- 
wald step rule. According to the conjecture of Stranski 
and TotomanowQ, the free energy barrier for the forma- 
tion of disordered fee should then be lower than the one 
for the formation of a CsCl cluster. The ST conjecture 
relies on the assumption that sub-critical nuclei are in 
quasi-equilibrium. However, as we will show, this turns 
out not to be the case. We can test this by repeating the 
FFS scheme with a different kinetic MC scheme that in- 
cludes an additional trial move: the swap of positive and 
negative particles. If the system is already in equilibrium, 
the introduction of additional MC moves will not change 
the structure of the sub-critical nuclei. 

When we performed FFS simulations including 20% 
swap moves we observe the formation of charge-ordered 
clusters with a CsCl structure (Fig[2]b). Not only the 
structure of the sub-critical nuclei has changed but also 
the size of the critical nucleus: it now contains 65 ± 15. 
Moreover, the probability that a solid cluster of 5 par- 
ticles will form a bulk crystal has increased to 10~ 15±1 . 
The fact that the pathway for crystal nucleation can be 
altered by artificially improving the sampling of config- 
urational space indicates that local equilibrium is not 



established during the natural nucleation dynamics (no 
swaps). This observation is in direct contradiction to the 
key assumption underlying the Stranski- Totomanow con- 
jecture. Our simulations suggest that the time it takes a 
cluster to grow from a small size to the critical size is too 
short to allow for efficient sampling of the accessible con- 
figurational space - as a result, it gets kinetically trapped 
in a metastable structure. We expect that such behav- 
ior will be common when there is no "easy" kinetic route 
(e.g. martensitic transformation) from the metastable to 
the stable crystal phase. 

At higher pressures (p* = 18), where the probability 
of forming a post-critical nucleus is higher and no FFS is 
needed to observe the transition, the same phenomenol- 
ogy is reproduced. In a kinetic MC simulation without 
swap moves a substitutionally disordered fee lattice is 
formed. The same is obtained in a Brownian Dynamics 
simulation. In contrast, when swap moves are included, 
the liquid transforms into a substitutionally ordered lat- 
tice. 
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FIG. 2: Radial distribution function of an 80 particles cluster 
obtained in a typical FFS path at T*=l and p*=15 compared 
with the radial distribution function of the solids CsCl and 
disordered-fec. (a) without swap moves, (b) with swap moves. 
The structure of the clusters varies with the way in which the 
configurational space is sampled. The unnormalised radial 
distribution functions of the clusters have been multiplied by 
20 in order to compare with the "bulk" radial distribution 
functions. 

Interestingly, as we can selectively prepare fee or CsCl 
nuclei by changing our kinetic MC scheme, we can now 
separately compute with umbrella sampling llj the free 
energies of these two different types of clusters. In 
umbrella-sampling simulations, the minimum free energy 
path is sampled along a given reaction coordinate. Nev- 
ertheless, when the calculation is carried out without 
swap moves, the growth of CsCl clusters is dramatically 
slowed down, instead, fee clusters are formed and per- 
sist for long time in the system. If the simulation is run 
long enough, the structure of the clusters changes into 
CsCl, suggesting that CsCl clusters have indeed lower 
free energy. By including swap moves in the umbrella 
sampling scheme, clusters grow directly in their lowest 
free-energy state (CsCl), yielding a different free-energy 
barrier. FigJ3] shows the barriers for both types of calcu- 
lations. As can be seen from the figure, the free-energy 
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barrier for the nucleation of fee clusters is higher than 
that for CsCl clusters. Snapshots of two typical critical 
clusters (fee and CsCl) are also shown in Fig. [3J 

This observation has direct implication for the inter- 
pretation of experiments 2^, 23|. The dynamics of real 



binary crystals of charged colloids is best described by 
the kinetic MC scheme (i.e. without unphysical swap 
moves) . Hence we should expect that in charged colloidal 
systems, crystallization proceeds through a sequence of 
non-equilibrium sub-critical crystal nuclei. In experi- 
ments on crystallization in binary charged colloids [l5| . 
both substitutional^ ordered (CuAu) and substitution- 
ally disordered (fee) crystallites have been observed. The 
simultaneous observation of both phases could be ex- 
plained thermodynamically if the experimental condi- 
tions fortuitously happened to correspond to coexistence. 
The present work suggests another explanation: non- 
equilibrium nucleation of the fee phase precedes a sub- 
sequent, slow transformation to the substitutional^ or- 
dered crystal phase. 
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tion inside the clusters is kinetically inhibited. A crystal 
cluster could change its internal structure by a succes- 
sion of particle additions and removals, but in practice 
this would mean that a disordered fee cluster would have 
to redissolve almost completely before it can form a CsCl 
cluster. The "success" of the small sub-critical fee clusters 
blocks the subsequent formation of the more stable CsCl 
clusters. This phenomenon is reminiscent of the "self- 
poisoning" of small crystallites during the rapid growth 
of post-critical crystal nuclei [24]. The difference is that, 



in the present case, the self-poisoning already takes place 
with sub-critical nuclei. 

The present results imply that, at least to predict crys- 
tal nucleation, there are situations where it is not enough 
to compute the free energy barrier that separates the par- 
ent phase from resultant solid structures - beyond a cer- 
tain cluster size, the formation of the lowest free-energy 
clusters may be kinetically inhibited. The fast growth of 
the clusters results in the breakdown of the local equilib- 
rium assumption for sub-critical nuclei. 
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organization for Scientific Research (NWO). 



% ' 50 100 ' 150 

n 



FIG. 3: (a) : Free energy barriers calculated with umbrella 
sampling at T*=l and p*=15 in a system of 8000 parti- 
cles. The CsCl clusters have lower free energy, but unless 
unphysical MC moves are used in the sampling, the system 
remains kinetically trapped in a disordered-fcc route of higher 
free energy, (b): Snapshot of a typical critical cluster with 
disordered-fcc structure and (c):CsCl structure. Note that 
the disordered critical cluster is bigger than the ordered one. 

At first sight, it might seem that the present results, 
although at odds with the ST conjecture, are not incom- 
patible with CNT. After all, within that theory, preferred 
nucleation of the crystal structure with the higher nucle- 
ation barrier is possible if a large kinetic prefactor in Eq.[T] 
compensates the effect of the higher nucleation barrier. 
Yet, the existing versions of CNT do not correctly de- 
scribe this effect: in CNT the kinetic prefactor describes 
the rate at which clusters grow due to the attachment 
and detachment of single particles to a pre-existing crys- 
tallite, and the rate of addition and removal of particles 
is hardly different for fee and CsCl clusters. What seems 
to happen is that small clusters have a disordered fee 
structure, but this structure cannot act as a template for 
subsequent CsCl growth, whilst a structural phase transi- 
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